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In hydrodynamical modeling of heavy-ion collisions, the initial-state spatial anisotropies are trans- 
lated into momentum anisotropies of the final-state particle distributions. Thus, understanding the 
origin of the initial-state anisotropies and their uncertainties is important before extracting specific 
QCD matter properties, such as viscosity, from the experimental data. In this work we review the 
wounded nucleon approach based on the Monte Carlo Glauber model, charting in particular the un- 
certainties arising from modeling of the nucleon-nucleon interactions between the colliding nucleon 
pairs and nucleon-nucleon correlations inside the colliding nuclei. We discuss the differences between 
the black disk model and a probabilistic profile function approach for the inelastic nucleon-nucleon 
interactions, and investigate the influence of initial-state correlations using state-of-the-art modeling 
of these. 

PACS numbers: 12.38.Mh,25.75.Cj,25.75.Ld 



I. INTRODUCTION 

In ultrarelativistic heavy ion collisions performed 
at Relativistic Heavy Ion Collider (RHIC) and Large 
Hadron Collider (LHC) elliptic flow - the second Fourier 
coefficient vi that quantifies the azimuthal anisotropy in 
the measured particle distributions - has been found to 
be large [l|-[|[. The appearance of a significant i>2 is a 
clear signature of pressure formation in the system. It 
is also consistent with predictions from relativistic hy- 
drodynamics @, IH : At non zero impact parameter the 
overlap area of the colliding nuclei is eccentric in the 
transverse coordinate plane. This spatial anisotropy is 
translated first into a flow- velocity anisotropy during the 
hydrodynamical evolution and, finally, at the decoupling 
of the system, into a measurable momcntum-anisotropy 
of final-state particle distributions. 

Recently it was found out that geometrical fluctuations 
in the nucleon transverse positions generate initial state 
anisotropies [f| also for odd harmonics v n . Nowadays, the 
experiments at RHIC and LHC have been able to mea- 
sure non zero flow coefficients up to vq ^ n recent 
years, there have also been considerable developments in 
event-by-event hydrodynamical modeling [Ill4l6| , thanks 
to which it has become possible to study the higher har- 
monics fl7rll9t and v\ at mid rapidity [20l - i22| . Studies 
based on hydrodynamics have shown that all initial-state 
anisotropies are transferred to the final measurable flow 
values in a similar way as eccentricity is translated to 
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elliptic flow. 

In the past few years much effort has been put on 
developing viscous hydrodynamics [23l - l36l ] and it has 
turned out that the elliptic flow and especially the 
higher harmonics arc sensitive to the (shear) viscosity 
r\ [16l ITfjl . l37l . [38| . Since the final flow observables closely 
reflect the initial-state spatial anisotropies, the uncer- 
tainties in the assumed (or computed) initial state are 
transferred to the computed final-state flow observables. 
In the studies where one has tried to estimate the shear 
viscosity-to-entropy ratio, 77/s, using elliptic flow mea- 
surements, one has seen that the initial eccentricity differ- 
ences between different initial state models lead to large 
uncertainties in the extracted value of 77/s [U IH[. For 
determining the shear viscosity from the flow measure- 
ments, it is very important to chart all the relevant un- 
certainties in the computation of the initial asymmetries. 

The Glauber model [39| is usually a key element in 
computing the initial states for hydrodynamical model- 
ing of ultrarelativistic heavy-ion collisions. Some years 
back, most hydrodynamical calculations assumed smooth 
initial states where the (energy or entropy) densities were 
assumed to scale with the density of binary collisions or 
wounded nucleons computed from the optical Glauber 
model; see, e.g. |40j . Now that the importance of the ini- 
tial density fluctuations has been realized, Monte Carlo 
Glauber (MCG) modeling has become more frequently 
used. So far the black disk (hard-sphere) modeling of the 
nucleon-nucleon (NN) interactions has been the stan- 
dard choice |4lTl43j in these studies, although also more 
involved probabilistic ways to model the NN interactions 
have been known for a long time [HI, I44l - l46j . 

For MCG modeling, one needs the nucleon configu- 
rations inside the colliding nuclei as an input for each 
event. In the simplest approximation the nucleons are as- 
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sumcd point-like and the nucleon positions are obtained 
by just randomly sampling the Woods-Saxon number 
density distribution. However, to improve the model- 
ing, TVTV correlations should also be considered, as sug- 
gested in Refs. (47l - [5l| . Experimental evidence for the 
two-body TV TV correlations generating high-momentum 
components of the nuclear wave function [52l . [53j . are 
discussed in Refs. [541 - 156} . Furthermore, recent studies 
in Ref. (5?| have shown that the effects from including 
more realistic TVTV correlations are not negligible e.g. in 
the generated initial eccentricity. 

In this paper we study for Au-Au collision at RHIC en- 
ergies, two different uncertainties in computing the initial 
state asymmetries from the MCG model: one related to 
the modeling of the inelastic TVTV collisions between nu- 
clcons from different nuclei, and one related to the TVTV 
correlations in the nucleon configurations in each of the 
colliding nuclei. The paper is organized as follows: In 
Sec. Ill Al we outline how the initial nucleon configurations 
can be obtained taking the TVTV correlations into account. 
After that, in Sec. Ill Bl we present two different ways to 
model the inelastic TVTV primary collisions. The studied 
initial-state anisotropies are defined in Sec. Ill CI and re- 
sults are presented in Sec. IIIII Finally, the conclusions 
are drawn in Sec. IIVI 



II. MONTE CARLO GLAUBER MODEL 
FRAMEWORK 

A. Nucleon configurations 

The initial state of a nucleus in the MCG calculations 
is usually taken as a collection of particles distributed 
according to a probability distribution given by the cor- 
responding (Woods-Saxon) number density distribution 
measured in electron scattering experiments. Given the 
complexity of the nuclear many-body problem, the effects 
of spatial, spin and isospin dependent correlations among 
the nucleons are usually overlooked and the nucleons are 
positioned randomly for each of the simulated events. 

Recently, in Ref. [47} it was shown how such an ap- 
proach can be modified by including initial states, which 
are prepared in advance, in the commonly used computer 
codes. Also, a method to produce such configurations 
was introduced. The method is based on the notion of a 
nuclear wave function which contains the nucleonic de- 
grees of freedom and which is used to iteratively modify 
the positions of randomly distributed nucleons using the 
Metropolis method so that the final positions correspond 
to the probability density given by |V>| 2 - The method 
is constructed to reproduce the same nucleon number- 
density distribution as the usual one and, in addition, to 
reproduce the basic features of the two-nucleon density 
in the presence of the TVTV correlations. The model wave 



function is taken in the form 

A 

ifj{x x , x A ) = Y[ fij <t>(xi, X A ) , (1) 

i<j 

where <j> is the uncorrelated wave function and fij are cor- 
relation operators; here, Xi denotes the position, spin and 
isospin projection of the i-th nucleon. The correlation op- 
erator contains a detailed spin-isospin dependence, which 
is as follows: 

fij = f C ( r ij) + f a ( r ij ) CT i ■ a j + f T ( r ij) T i ■ T ] 

+ r {Tufa ■ ■ Tj + f (n 3 )S %J (2) 

+ f tT {r ij )S ij r i -T j + ..., 

where ry is the relative distance between nucleons i and 
j, cri and Ti are the Pauli spin and isospin operators, 
respectivel y, a nd Sij is the so-called tensor operator (see 
e.g. Ref. [52j ) which depends on the spin and spatial 
variables of nucleons i and j. In the most general case, 
the sum extends over a number of channels that are the 
same as the one appearing in modern nucleon-nucleon 
potentials used to successfully describe a variety of prop- 
erties of light and medium-heavy nuclei within different 
ab initio approaches. When all the correlation functions 
in Eq. ©, except the first one [f c {r)], vanish, we have 
the central correlation case. 

In this paper we will consider configurations produced 
both with the central correlations only (denoted in the 
following by central) and with a full set of the six correla- 
tion functions of Eq. © (denoted by full), including the 
tensor operator. These correlation functions }52j were 
developed by variational method for nuclei lighter than 
197 Au considered in this paper; nevertheless, they repre- 
sent the best approximation that we can offer, since no 
corresponding calculations exist for heavy nuclei. More 
specifically, we used correlation functions obtained for 
40 Ca. Many recent theoretical works support the univer- 
sality of TVTV correlations in nuclei |58l - l60l ] . We justify the 
use of 40 Ca correlation functions by the observation that 
they differ very little from those obtained for 16 O; 40 Ca 
is a large enough nucleus to neglect additional differences 
with 197 Au due to A dependence in this context, as well 
as differences arising from the fact that correlation func- 
tions were obtained for doubly magic, symmetric nuclei 
such as 16 and 40 Ca. Moreover, as shown in Ref. [48| . 
the two-body densities resulting from configurations ob- 
tained using these correlations are clearly more realistic 
than the completely uncorrelated ones. 

The method developed in Ref. (4?J allows in princi- 
ple to deal with any set of correlations; nevertheless, 
although the state dependent correlations introduce an 
extra computational effort due to their non commuta- 
tive nature, they can be treated up to a certain degree. 
In this paper we have used configurations including two- 
body full correlations, and configurations of three-body 
clusters surrounding each of the nucleons, induced by full 
correlations. Genuine three-body correlations will also be 
discussed in the following. 
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B. Modeling the inelastic interactions 



Thus, we have 



We work in the Glauber model framework [39J , neglect- 
ing the effects of inelastic diffraction that lead to fluctua- 
tions of the strength of the TV TV interactions [6l|, HH . To 
generate the inelastic TV TV collisions of interest here, we 
use the following two different approximations for decid- 
ing whether a collision between the nucleons i and j from 
different nuclei takes place: 

• Black disk approximation, used recently, e.g., in Rcf. 
[l5| . where one assumes the two nucleons to interact in- 
elastically with a probability one if their transverse sepa- 
ration bij is within a radius defined by the inelastic TV TV 



a 



NN 



(3) 



• Profile function approach, where the probability of an 
inelastic interaction between the nucleons i and j is given 
by 



p(6«) = i - |i-r(&y)r 



(4) 



and where the profile function T is expressed in terms of 
the total and clastic TV TV cross sections as follows: 



r(b, 



T NN -bl/(2B) 



4ttB 



(5) 



with B = {a^ N f/{lQ-Ka NN ). 

The probability distribution P(bij) can be derived in 
the Born approximation of the potential-scattering for- 
malism, by parametrizing the TV TV elastic scattering am- 
plitude as [Hi, m, E3 



= C(i + a)k e _, Bq , 
Air 



(0) 



where C and B are constants to be determined, k ~ E 
and a as for ultra-rclativistic energies. For the small- 
angle scatterings of interest here, the vector q corre- 
sponds to the difference between the incoming and scat- 
tered wave vectors in the transverse plane, and q = |q| 
relates to the scattering angle as q as kO. We fix C to the 



measured er^Jv 011 tne basis of the optical theorem, 
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*Tn = T lm[f(0)]- 



(7) 



The profile function in Eq. (|5|) is obtained as a Fourier 
transform, 
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and the elastic TV TV cross section from this as 
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(10) 

(11) 



and B can be fixed on the basis of measured cross sec- 
tions. For the current set-up at ^Snn = 200 GeV, we 
take B = 14 GeV~ 2 and er^v = 52 nib, which correspond 
to cr^jy = 42 mb and o~^ N = 9.9 mb. 

Finally, using Eqs. (|7|), © and (fTT|) . we arrive at the 
probability function of Eq. ((4]), whose integral over the 
transverse separation gives the inelastic TV TV cross sec- 
tion: 



'NN 



d 2 b l3 (2ReT(b l3 ) - |r(6 tf )| : 



d'bu i - |i - r(6 tf )l' 



(12) 



C. Spatial asymmetries and their fluctuations 

In this section we define the initial state anisotropics, 
which are studied in this work. We focus on the first 
three harmonics n = 1,2,3: dipole asymmetry, eccen- 
tricity, and triangularity. Higher harmonics are left out 
since they are more complicated due to the fact that they 
can be mainly originating from the lower harmonics [63J . 
Also, from the experimental results we know that the 
second and third harmonics are the largest ones. 

We calculate the asymmetries from the wounded nu- 
cleon positions which are obtained from the MCG model. 
In the following, the angle brackets denote an average 
over wounded, or participant, nucleons. The asymme- 
tries arc defined as 



(w(r) cos(n(</> — ipn)} 
(w(r)) 



(13) 



where w(r) is a weight and ip n is an orientation angle 
that is obtained as 

1 (w(r) sm(nd))) ir 
^-arctanf-^ y%L + -, (14) 

n {w(r) cos(n<p)) n 

where arctan is always placed in the correct quadrant. 
These quantities are always calculated in the center-of- 
mass of the wounded nucleon system. We choose to use 
w(r) = r 3 for n = 1, w(r) = r 2 for n = 2, and w(r) = r 3 
for n = 3 0. 

We also study the fluctuations of the initial-state asym- 
metries since the current flow analysis methods are sen- 
sitive to the fluctuations of the flow coefficients 641. 



Since final flow values reflect the initial-state asymme- 
tries, the fluctuations of the flow coefficients should follow 
the initial-state fluctuations. We define the fluctuations 
of the anisotropies as 



Ae r 



TV 



(15) 
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where t % n is the asymmetry in event i and (e„) denotes 
the average over TV events. 



III. RESULTS 

A. Number of wounded nucleons and binary 
collisions 

First we look how the interaction models and TV TV cor- 
relations affect the number of the participants (-/V par t) 
and binary collisions (N co \\). In Fig. [1] we have plotted 
these as a function of the impact parameter. One can 
see that in central collisions the results are the same with 
both interaction models, but when impact parameter be- 
comes large, b = 10 — 15 fm, i.e., when only the edges of 
the nuclei collide and when the number of participants 
and binary collisions are of the same order, the differ- 
ence between models is over 5 %. This holds both for 
wounded nucleons and binary collisions. From Fig. [T] it 
can be seen that the effects of TV TV correlations on these 
quantities are very small. 

B. Effects of the interaction model on e n and their 
fluctuations 

Next, we consider the two different interaction mod- 
els discussed in Sec. Ill Bl and present their effects on the 
anisotropies defined above. In these calculations the TV TV 
correlations are neglected for clarity; the effects of these 
correlations will be studied in the next section. In Fig. [2] 
we have plotted the anisotropies £i, £2, and £3 as a func- 
tion of iVp ar t using the black disk and profile function ap- 
proaches to model the inelastic TV TV collisions. The black 
disk approximation results in larger dipolc asymmetry ei, 
eccentricity £2, and triangularity £3. In all cases, the ob- 
tained asymmetries in the profile function approach are 
slightly (< 10 %) smaller than in those in the black disk 
case. Apparently, toward peripheral collisions, we have 
more fluctuations and a less-well-defined shape with the 
profile function approach, hence a smaller £„. 

The difference between these two cases is negligible in 
most central collisions for every n. When moving toward 
peripheral collisions, the results start to deviate from 
each other and the largest difference is approximately 
on the order of 10%. In central collisions most of the nu- 
cleons experience several collisions and, thus, the details 
of the collision model are not very important. In periph- 
eral collisions, however, more and more of the nucleons 
collide only once or twice, meaning that the interaction 
model details start to play a role. One must also remem- 
ber that as we saw from Fig. [TJ the events with the same 
impact parameter are mapped to slightly smaller TV part 
values with the black disk interaction model. 

In Fig. [3] we plot the relative fluctuations of the initial- 
state anisotropies. Here the order of the curves is oppo- 
site compared to the e n when n = 2, 3. This means that 
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FIG. 1. (Color online) Number of wounded nucleons (a) and 
binary collisions (b) as a function of impact parameter with 
different interaction models and initial-state configurations. 
The results shown correspond to the black disk approxima- 
tion with uncorrelated configurations (green dashed line) ; the 
profile function approach T(bij) with uncorrelated configura- 
tions (black solid line); the profile function approach with 
central NN correlations (red squares); and the profile func- 
tion approach with full NN configurations (blue circles). The 
insets show the ratios of the last three cases to the black disk 
+ uncorrelated one, with corresponding notations. 



at least partly the difference in the relative fluctuations 
is explained simply by the fact that the absolute value 
is slightly higher in the case where the relative fluctu- 
ations are smaller. For e± the relative fluctuations are 
larger with black disk interaction, indicating that then 
also the absolute fluctuations are larger with black disk 
than profile interaction. 

The relative fluctuations of £1 show a decreasing trend 
when approaching central collisions but the impact- 
parameter dependence is relatively weak. On the other 
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FIG. 2. (Color online) First three harmonics e„ as a function 
of the number of participants. Results shown are for the two 
different interaction models. In these plots no NN correla- 
tions are taken into account. 



FIG. 3. (Color online) Relative fluctuations of the first three 
harmonics e n as a function of number of participants. Results 
are shown with two different interaction models. In these 
plots no NN correlations are taken into account. 



hand, the relative fluctuations of €2 clearly depend on 
the impact parameter. In central collisions the fluctu- 
ations are large since the system is azimuthally very 
symmetric, and toward peripheral collisions the fluctu- 
ations decrease, since the collision area becomes clearly 
eccentric. Triangularity £3 has yet another behavior: its 
relative fluctuations stay approximately constant when 
Npart > 100. Since the triangularity is created purely by 
the fluctuations in the positions of wounded nucleons, we 
can expect that the relative fluctuations have no central- 
ity dependence. All in all, Aei/ei and Ae^/es have only a 
mild centrality dependence, while A^j exhibits a clear 
dependence on the centrality. In addition, we conclude 
that for these relative fluctuations, the NN interaction 
model uncertainties do not play a major role. 



C. Effects of the NN correlations on e n and their 
fluctuations 

Next, we investigate the effects of different models for 
the initial state NN correlations. We have chosen to use 
the r(&jj) collision profile function in these calculations. 
In Fig. 2] we have plotted the obtained anisotropics with 
three different NN correlations: no correlations at all, 
only central correlations and full correlations. From all 
panels we see that no correlations and full correlations 
are very close to each other, but the central correlations 
have a smaller anisotropy. 

In all cases the difference between central correlations 
and the two other cases is largest at central collisions 
and it gets smaller toward peripheral collision. Nuclcon 
correlations are most important in the middle of the nu- 
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FIG. 4. (Color online) First three harmonics e n as a function 
of number of participants. Results are shown with different 
NN correlations and using the T(bij) collision profile function 
approach. 



FIG. 5. (Color online) Relative fluctuations of the first three 
harmonics e n as a function of the number of participants. 
Results are shown with different AW correlations and using 
the T(bij) collision profile approach. 



cleus since the nucleon density is largest there. Thus, 
in the central collisions, where most of the wounded nu- 
cleons come from the middle of the nuclei, the effect on 
anisotropics is largest. When moving toward the periph- 
eral collisions, the relative amount of wounded nucleons 
coming from the edges of nuclei is increasing. Thus, in 
very peripheral collisions where only the edges are over- 
lapping, the effect disappears. 

We must bear in mind that central correlations have 
only repulsive character; their effect can be mimicked 
by an ad-hoc exclusion radius, as shown in Ref. [57| . 
Full state-dependent correlations, instead, have a com- 
plex structure, which causes nucleons to be at a given 
average distance in the nucleus and results in nuclear 
binding and its saturation with increasing A. Partially 
including these full correlations and disregarding three- 
body repulsion, as it is explained below, produces the net 
result of working in the opposite direction of repulsion, 



for the considered quantities. 

Next, in Fig. [SJ we plot the relative fluctuations of 
anisotropics with different initial state NN correlations. 
We can see that inclusion of full correlations brings the 
results back toward the uncorrelated case; this is less ev- 
ident for Ae3/e3. However, now the difference is small in 
the central collisions and it is largest in the semi periph- 
eral collisions. The difference again vanishes at the most 
peripheral collisions. 

The results obtained with realistic configurations de- 
serve some discussion. The production of configu- 
rations with full realistic correlation functions differs 
substantially from the central correlations case, since 
in the realistic description there are several spatially- 
dependent correlation functions complemented with spin- 
and isospin-dependent operators, as shown in Eq. 
At this stage, we have performed a truncation of the 
chains induced by realistic correlations at the level of 
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three-particle chains; in principle, due to the non com- 
mutativity of the operators in Eq. ([2)1. the chains are 
A-body operators. Moreover, we restricted our calcu- 
lations to three particles which are within a given ra- 
dius from the active particle in the Metropolis algorithm. 
It should be stressed that these kinds of simplifications 
are not suggested by specific physics arguments but are 
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FIG. 6. (Color online) The effect of different models of realis- 
tic correlations on the relative fluctations of eccentricity; see 
text for explanation. 



rather dictated by the enormously increasing computing 
time. As a result, the truncation obviously induces some 
amount of uncertainty in our results. The overall trend, 
that the full correlations seem to bring the results back 
toward the no correlation case, is nevertheless what we 
wish to emphasize here. 

In order to illustrate the effect of the truncation, we 
have repeated the calculation with configurations includ- 
ing only two-body, realistically correlated clusters (2b 
only) , and compared with the three-body calculation out- 
lined above (3b chains). The results for Ae2/e2 are com- 
pared with the corresponding uncorrelated and central 
correlation case in Fig. [5] Interestingly, it can be seen 
that the full correlations with the 2 body chains cause 
an effect into the opposite direction than the more ad- 
vanced 3b chains do. Although beyond the scope of this 
paper, it would be interesting to study how sensitive the 
considered anisotropies, and also the higher moments of 
participant matter distribution, are to the higher-order 
chains of two-body state-dependent correlations and to 
the genuine three-body correlations which we have not 
included here. The outcome of such a study is, however, 
difficult to predict, due to the complicated interplay of 
attraction and repulsion between the three particles in 
different spin and isospin states. 



IV. CONCLUSIONS 

In this paper we have charted some of the uncertainties 
in the computation of the initial state anisotropies from 
the Monte Carlo Glauber model. We used two differ- 
ent ways of modeling the inelastic interactions between 
the colliding nucleons. The difference between these two 
cases gives us an estimate about the uncertainties related 
to this part of the model: in central collisions the details 
of the interaction model play a minor role, but in the 
peripheral collisions such details can cause uncertainties 
up to 10% in the first three harmonics e\, 62 and £3. We 
also checked that with these two interaction models the 
difference in the number of wounded nucleons and binary 
collisions remains small in central collisions, but at im- 
pact parameters 10-15 fm the difference can be around 
10%. We also note that during the writing process of this 
article, a similar nucleon interaction model study was re- 
leased in Ref . [|35[ . The main differences to our study are 
the different form of the elastic NN scattering amplitude 
as well as the treatment of the NN correlations. 

We also presented a study of the effects of NN corre- 
lations with an update of correlated configurations and 
extended discussion as compared with the previous pub- 
lished papers on this subject. We confirmed that the 
inclusion of centrally correlated nucleon configurations 
produce the effects to eccentricity and its relative vari- 
ance as was claimed by Ref. [13] ■ As a new result, we ob- 
served that the inclusion of realistically correlated config- 
urations (two-body full correlations, three-body chains) 
seems to essentially cancel this effect and bring the re- 
sults back close to the no correlations case. The effect 
is similar for dipole asymmetry and triangularity as for 
eccentricity. However, we also showed that there are still 
uncertainties caused by the truncation done in the nu- 
cleon configuration calculation with full correlations and 
we expect three-body correlations to play a role. 

In this paper we studied two sources of additional un- 
certainty in the Monte Carlo model calculations for the 
initial state anisotropics. The uncertainty caused by the 
studied effects to these anisotropics was found to be max- 
imally of the order of 10%. Now that - thanks to the 
recent developments in event-by-event hydrodynamics - 
more precise comparisons of flow coefficients between the 
data and the theory are becoming possible, it is impor- 
tant to chart all the relevant uncertainties to this preci- 
sion, so that the QCD matter properties could eventually 
be determined from the measured particle spectra and 
their azimuthal asymmetries. 
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